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Abstract The green odorous frog (Odorrana 
margaretae) has an interesting ring-shaped divergence 
pattern around the Sichuan Basin and documenting its 
morphological variations is essential in understanding 
its evolutionary history. Using curvilinear models, 
we detected significant geographical clinal variations 
in morphological traits, particularly sizes, of female 
O. margaretae. Males had significantly smaller sizes 
than females, and also had smaller variation ranges 
than females. One major trend of morphological 
variations was clinal: populations from the west 
tended to have a larger size with wider head and 
longer posterior limbs than populations from the east. 
Species history, with an early extended isolation and 
two subsequent secondary contacts, may explain most 
of the geographical clinal variations of O. margaretae. 
Bioclimatic factors may also contribute in explaining 
the variance of morphology. 
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1. Introduction 


The odorous green frog (Odorrana margaretae) displays an 
interesting ring-shaped divergence pattern around the Sichuan 
Basin of western China, much like a ring species (Qiao et al. 
2018). It is a large stream-dweller primarily distributed in the 
mountains of western China with a few sporadic distribution 
records at the east (Figure 1A; Fei et al., 2009). Using DNA 
sequence and microsatellite DNA data, Qiao et al. (2018) 
examined its phylogeographical history. The current ring- 
shaped distribution pattern likely originated from two refugial 
populations, one at the west and the other at the southeast of 
the Sichuan Basin. Both populations expanded around the Basin 
and formed two contact zones. Extensive admixture occurred 
at the south contact zone, which became an evolutionary 
‘melting pot’, and the second contact zone at the northwestern 
corner of the Basin only revealed limited hybridization and 
partial reproductive isolation has developed between the two 
expansion fronts (Qiao et al., 2018). Furthermore, the chain 
populations demonstrated a strong isolation by distance 
pattern around the Basin, suggesting that the genetic variation 
were mostly gradual and continuous (Qiao et al. 2018). 
Morphological variations among these populations have 
been noticed (Fei et al, 2009; Shen et al., 2009). For example, 
Fei et al. (2009) reported that individuals from the western 
populations (Mt. O'mei) are larger than individuals from the 
eastern (Wushan) and southern (Anlong) populations. Large 
ventral color pattern variations were also described (Fei et 
al. 2009). With its ring-shaped divergence pattern, we would 
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Figure 1 A: Sampling sites of Odorrana margaretae around Sichuan Basin. B: Patterns of geographical clinal variation in male and 
female morphology: mean scores of PC1 from every site were plotted in coordinates. The grey scale and size of symbols are according 


to increasing size. 


expect morphological clinal variations around the Sichuan 
Basin, in a similar fashion to genetic variation. Numerous early 
works have demonstrated that both primary differentiation 
and secondary contact can result in clinal variation (eg. Endler 
1977), and clinal variation in morphological traits has been 
well documented in various animal groups (e.g, mammals, 
Okeefe et al, 2013; John and Richmond, 1993; turtles, Ennen et 
dl., 2014; geckos, Fitness et al., 2012; frogs, Poynton and Loader, 
2008; fishes, Moore and Hendry, 2005). Understanding the 
significance and the mechanisms of clinal variation will lead to 


a broad understanding of fundamental evolutionary processes 


such as local adaptation, population differentiation, and 
ultimately, speciation (Mayr, 1942; Endler, 1977). 

In this study, we examined the morphological traits of 
O. margaretae populations around the Sichuan Basin. Our 
objectives are 1) documenting morphological variations 
quantitively and 2) exploring potential causes of these 
variations. We specifically considered two alternative causes, 
its evolutionary history in particular the formation of its ring 
distribution, or climatic factors. If the evolutionary history is 
the leading cause, we would expect a clinal variation follow 


the ring distribution; if climate dominates these variations, we 
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would expect a strong correlation between morphological traits 


and climatic factors, irrelevant to its ring distribution. 
2. Materials and Methods 


21. Data collection A total of 152 adult O. margaretae specimens 
(61 females, 91 males) from 20 locations around the Sichuan 
Basin were measured (Appendix Table SI, Figure 1A). Specimens 
were collected from 1956 to 2018 (Appendix Table S1), and 
are deposited at the Herpetology Collection of the Chengdu 
Institute of Biology, Chinese Academy of Sciences (in 1096 
formalin), and at the Henan Normal University. We selected 
21 external body characteristics that are most commonly used 
for anurans (Table 1). All measurements were collected using 
a Mitutoyo ditgimatic caliper (Mitutoyo Corp., Japan) and 
measured to the nearest 001mm. 

For environmental factors, 19 bioclimatic variables for the 
years 1970-2000 with a spatial resolution of 30 arc-seconds 
resolution (-1km) were first obtained from WorldClim version 
2 database (http://www.worldclim.org/; Fick and Hijmans 
2017), and extracted by site using ArcMap 102. Eleven of these 
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variables are related to variations in temperature (BIO1-11, 
Table 2), and eight are related to variations in precipitation 
(BIO12-19, Table 3). 


2.2. Data Analysis We first compared the snout-vent length 
(SVL) between females and males using ANOVA. Since 
significant differences were detected, all downstream analysis 
was conducted separately for males and females. 

The 21 morphological measurements were first subjected 
to a principal component analysis (PCA). For morphological 
traits, we plotted the mean scores value of PCI from every 
site. Wilcoxon rank sum test was used to compare the 
morphological data among different geological groups by using 
the first principle component score values. The Wilcoxon tests 
were conducted on an online calculator (https://astatsa.com/ 
WilcoxonTest/). 

To avoid autocorrelation, a PCA was also performed on 
the 19 bioclimatic variables. Values of bioclimatic variables 
with loading score greater than 02 on PCI and PC2 were then 
plotted around the Basin. The PCA analyses were carried out 
with the 'prcomp' function using R (v. 3.6.1), and plotting was 
conducted using R package ‘ggplot2’ (Wickham 2016). 


Table 1 Summary of 21 morphological characteristics and the corresponding PC1 (PC1morph) loadings on PCA in Odorrana margaretae. 
All measurements are in unit of millimeter. Values in bold represent loading absolute values greater than 0.2. Asterisks indicate significant 


difference between the sexes. 


Body measure Ae = 
Mean (+SD) PC1morph (79.37%) Mean (+SD) PC1morph (90.57%) 

Snout-vent length, SVL* 72.75 (£5.24) 0.6205 87.67 (£10.21) 0.6586 
Head length, HL 22.04 (+1.98) 0.1817 25.83 (£2.35) 0.1236 
Head width, HW 24.39 (£1.97) 0.2194 30.35 (+3.91) 0.252 
Snout length, SL 9.44 (+0.88) 0.0541 11.23 (+1.22) 0.0546 
Internarial distance, IND 8.61 (+0.68) 0.059 10.37 (+1.09) 0.0596 
Interoptic distance, IOD 5.93 (+0.75) 0.0388 6.87 (+1.26) 0.0563 
Antoptic distance, AOD 12.66 (+1.03) 0.1041 15.26 (+ 1.88) OBS 
Postoptic distance, POD 19.28 (+1.52) 0.155 22.72 (+2.54) 0.1536 
Eye-nostril distance, EN 5.20 (+0.45) 0.0363 6.11 (£0.75) 0.0353 
Eye diameter, ED 8.72 (£0.71) 0.0504 9.82 (+0.87) 0.0404 
Tympanum diameter, TD 3.91 (+0.37) 0.0178 4.28 (+0.46) 0.0154 
Hand length, HAL 18.29 (+1.50) 0.1457 21.50 (+2.51) 0.138 
Forearm length, FLL 21819. (3:1:5])) 0.1589 25.04 (+2.68) 0.1635 
Forearm width, FLW 7.71 (£1.27) 0.0691 8.14 (+1.26) 0.0352 
Foot length, FL 42.93 (£3.36) 0.3688 50.99 (£5.32) 0.3335 
Tarsus length, TSL 20.68 (+1.56) 0.1522 24.30 (+2.58) 0.1569 
Inner metatarsal tubercle length, IMTL 4.28 (+0.51) 0.0431 5.15 (£0.74) 0.0378 
Inner metatarsal tubercle width, IMTW 1.96 (+0.28) 0.0138 2.32 (+0.37) 0.0148 
Tibia length, TL 44.52 (+2.99) 0.3476 52.56 (+5.32) 0.3414 
Tibia width, TW 10.58 (+1.25) 0.1136 13.17 (+2.08) 0.1223 
Thigh length, THL 42.35 (£3.28) 0.3752 50.03 (+5.50) 0.3494 
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Table 2 The loadings on the first two components (PC1bioc and PC2bioc) of PCA of 19 bioclimatic factors. Values in bold represent loading 
absolute values greater than 0.2 (the cutoff was inferred from Sidlauskas et al., 2010). 


Bioclimatic factors PC1bioc (86.36%) PC2bioc (12.48%) 
BIO1 _ Annual Mean Temperature -0.0018 -0.0066 
BIO2 _ Mean Diurnal Range (Mean of monthly (max temp - min temp)) -0.0021 -0.0104 
BIO3 _ Isothermality (BIO2/BIO7) (* 100) -0.0069 -0.0376 
BIO4 _ Temperature Seasonality (standard deviation *100) 0.0641 0.3896 
BIO5 _ Max Temperature of Warmest Month -0.0019 -0.0047 
BIO6 _ Min Temperature of Coldest Month -0.0015 -0.0058 
BIO7 _ Temperature Annual Range (BIO5-BIO6) -0.0004 0.0011 
BIO8 _ Mean Temperature of Wettest Quarter -0.0012 -0.0025 
BIO9 _ Mean Temperature of Driest Quarter -0.0023 -0.0105 
BIO10 _ Mean Temperature of Warmest Quarter -0.0007 -0.0007 
BIO11 .. Mean Temperature of Coldest Quarter -0.0023 -0.0105 
BIO12 _ Annual Precipitation -0.7526 0.5832 
BIO13 _ Precipitation of Wettest Month -0.1891 -0.2453 
BIO14 _ Precipitation of Driest Month -0.0165 0.0459 
BIO15 _ Precipitation Seasonality (Coefficient of Variation) -0.0032 -0.1301 
BIO16 _ Precipitation of Wettest Quarter -0.4309 -0.4373 
BIO17 _ Precipitation of Driest Quarter -0.0543 0.1523 
BIO18 _ Precipitation of Warmest Quarter -0.4492 -0.4354 
BIO19 _ Precipitation of Coldest Quarter -0.0543 0.1523 


Table 3 Regression analysis for morphology, geography and bioclimatic factors. Six regression models are compared. (1) The cubic 
polynomial regression model that morphology (PC1morph) is regressed onto geographical coordinates (x, longitude; y, latitude); (2) (3) The 
linear regression model that morphology is regressed onto bioclimatic factors (PC1bioc and PC2bioc); (4) The linear regression model that 
morphology is regressed onto locality altitude; (5) The linear regression model that morphology is regressed onto the year of collection; (6) 
The full regression model that morphology is regressed onto geographical coordinates as well as bioclimatic factors. Significant predictors are 
indicated by “*” (P < 0.05). R^ is the correlation coefficient between the outcomes and the predictors, and P-value represents significance of the 
regression analysis. Akaike information criterion, AIC, was used to compare the goodness of fit of the models. 


Male Female 
id n R(X) P-value AIC oe R(%) P-value ^ AIC 
UPN E 0.8393 02167 8132 E 07725 00414 14375 
(2) PCImorph-PCIbioc : 0.0021 0.8758 90.88 : 0.0008 0.91 155.87 
(3) PCImorph=PC2bioc PC2bioc* ^ 0.6166 0.0000 7749  PC2bioc* 0.3397 0,0088 148 
(4) PC1morph-alt : 0.1271 0.2108 8901 - 00199 0.5651 155.5 
netas E 00133 0.6951 9073 s 0.0389 04187 15513 
(6) PC1morph-polym (x,y) + PC2bioc — 0.8402 0.3901 83.24 y" xy" xy" 08224 0.038 141.05 


For testing correlation between morphological traits and 
geographical and climatic factors, six models were analyzed and 
compared (Table 3). First, a trend surface analysis (Borcard et 
al, 1992; Legendre and Legendre, 1998; Botes et al., 2006; Cardini 
et al., 2007) was applied to fit geographical coordinates to 
variations in morphology, taking into account of nonlinearities. 
The morphological variables were regressed onto a third-order 


polynomial of longitude and latitude: 


Fx, y) = bo + bix + bzy + b3x? + baxy + bsy? + box? + byx?y + baxy? + boy, 
where x and y are longitude and latitude respectively. Second, 
a basic ordinary least squares (OLS) model was used for 
fitting bioclimatic factors, locality altitude, and the year of 
collection to morphological variations separately. Finally, a full 
model evaluation, including geographical variables and other 
variables that were identified as significant predictors in the 


OLS model, was conducted to explain variance in morphology. 
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The first principle component of PCA on morphology and the 
first two principle components of PCA on bioclimatic factors 
were included in relevant models. Akaike information criterion 
(AIC) was used to compare the goodness of fit of the models. 
All correlation models were conducted separately on males 
and females. The ‘Im’ function in R was used to archive all the 


regressions. 
3. Results 


Morphological data revealed significant sexual size dimorphism 
in O. margaretae. The SVLs of females were significantly larger 
than males (Table 1). The raw measurement data are provided 
in Appendix Table S2. 

The principle component analysis on 21 morphological 
measurements revealed that the first axis explained the vast 
majority of total variations in both males (PCI, 79.37%; PC2, 
5.28%) and females (PCI, 90.57%; PC2, 3.33%). For both males 
and females, PCI represented mostly body length (SVL), head 
width (HW) and hind leg length (FL, TL, and THL) (loading » 
02; Table 1). The highest loading was for snot-vent length (SVL) 
(males: 0.6205; females: 0.6586). 

The males and females of O. margaretae displayed similar 
variation patterns of morphology around the Sichuan Basin, 
but the pattern was stronger in females than in males (Figure 
1B). Individuals at the eastern region of the Basin (site 7-8, 
Figure 1) had the smallest body size, with mean PCI scores of 
-12.03 (males) and -26.37 (females). The body size increased both 
westward and northward around the Basin, but the western 
populations (site 14-20) attained a much larger size (mean PCI 
scores: males, 4.28; females, 9.13) than the northern populations 
did (site 1-6, the mean PC scores: males, -0.98; females, -11.04), 
with a significant difference between those two regions (males, 
Wilcoxon W = 0.99, P < 0.05; females, Wilcoxon W = 3.06, P < 
0.005 ). The southern populations possessed an intermediate 
size (site 9-13, the mean PCI scores: males, 0.99; females: 3.06) 
between the eastern and western populations (Appendix Table 
S3). No significant difference was found between the southern 
populations and others. 

PCA of 19 bioclimatic factors summarized 98.84% of the 
variation in the first two components, with PCI explaining 
86.36% of total variance, and PC2 explaining 12.48% (Table 2). 
PCI mostly represented the precipitation related factors (e.g, 
BIOD, BIO16 and BIO18), with the highest loading for annual 
precipitation (BIOI2). PC2 was correlated positively with annual 
precipitation (BIO12) and temperature seasonality (BIO4), and 
moreover, was negatively correlated with several precipitation 
factors (e.g., BIOI3, BIO16 and BIO18). For precipitation, the 
annual precipitation was high in both eastern and western 
regions, and became drier toward the north. The west of the 
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basin showed noticeably abundant rainfall in the wettest 
month/quarter and warmest quarter (Appendix Figure S1A). 
As for temperature, the western region possessed a more stable 
seasonality than others; the temperature seasonality reduced 
gradually from east to west (Appendix Figure SIB). 

We detected a significant correlation between geography 
and morphology in females, and between bioclimatic factors 
and morphology in both sexes (Table 3). The best-performing 
model in females was the full regression model that included 
both geographical and bioclimatic factors (with the minimum 
AIC: 141.05), and the spatial components (y, xy and x'y) were 
significant in explaining the morphological variations. The 
model that included only geographical components also 
provided good fit for morphology (AIC: 14375), and geography 
explained 77.2596 the variance of the morphological variations 
in females (Table 3). Male morphology was not significantly 
correlated with geography, but was significantly correlated 
with PC2 of the bioclimatic factors (F; = 9.547, P < 0.005, Table 
3). Neither locality altitudes nor collecting years showed any 
significant effects on morphology. 


4. Discussion 


There are substantial morphological variations within the 
green odorous frog around the Sichuan Basin. The variations 
exhibit a clear geographical clinal pattern. While the eastern 
populations possess the smallest body size, the size increases 
along the southern margin of the Basin and the western 
populations attain the largest body size. Along the northern 
margin, populations also increase their body sizes, albeit to a 
less degree compared to the western populations (Figure IB, 
Table 1). Furthermore, sexual size dimorphism is significant; 
females are larger than males and also show a larger variation 
range than males. For example, females have an SVL range of 
6796-109.69 mm, while the range is 58.89-84.67 mm in males 
(Appendix Table S2, Table 1). 

Climate has a clear effect on morphological traits of 
this species. We detected a significant correlation between 
the morphological traits and PC2 of climate data (Table 3). 
Nevertheless, the climate data did not act as a significant 
predictor in the full model; this is likely because of the 
commonly observed autocorrelation between climate 
and geography. Climate often has significant influence on 
phenotypic variations both plastically or genetically, which 
is well-documented (e.g. Hubbe et al., 2009; Siepielski et al., 
2017; Urban et al., 2014). The western region of the Sichuan 
Basin has a warm and stable climate that is wet overall and 
receives abundant rainfall during the wettest period of the year 
(Appendix Figure SI). These favorable environmental factors 
likely promote growth. Furthermore, high environmental 
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humidity, long wet periods, and mild winters often improve 
larval survivorship and breeding success of those aquatic 
breeders (Banks et al. 1994; Blaustein et al, 2010; Scherer et al, 
2008). All of these may have contributed to the large body size 
of green odorous frogs in this region. The initial isolation of 
the species (Qiao et al. 2018) likely produced a west-large and 
east-small dichotomy, and during the westward expansion 
of the east refugial population along both the southern and 
the northern margins of the Basin, the body sizes increased. 
Although this parallelism could be caused by several other 
mechanisms (eg. chance, boundary effect), climatic factors likely 
played an important role. The observed climatic variations 
may have also caused several other phenotypic variations of 
the species. As important cues for the onset of reproduction, 
high temperature and precipitation may also promote an 
early breeding season, which is likely beneficial (Altwegg and 
Reyer, 2003). The breeding season of the western populations 
are generally earlier than the eastern populations. For example, 
the Mt. O'mei population at the west starts breeding during 
winter or early spring (Fei et al, 2009), while the Mt. Tian Ping 
population at the east has a breeding season from late August 
to early September (Shen et al. 2016). 

Species evolutionary history may also explain the observed 
morphological variations. The observed clinal variation is 
consistent with its ring distribution pattern and is compatible 
with the two-refugia hypothesis of its species history (Qiao 
et al. 2018). The green odorous frog likely evolved from two 
historical refugia formed during Pleistocene, the western region 
and the south-eastern region of the Sichuan Basin (Qiao et al., 
2018). The extended isolation period, which produced the initial 
genetic divergence, may have also generated the morphological 
difference, mostly body size difference with the eastern 
population being small and the western population being large. 
Climatic differences induced natural selection likely contributed 
to the observed phenotypical variations, while genetic drift 
likely also played an important role, as it has been demonstrated 
in many cases (eg, cichlids, Arnegard et al., 1999; Van Oppen 
et al., 1997; Markert et al., 1999). Post-glacial expansion and 
subsequent hybridization likely produced the geographical 
clinal variation (Figure IB). The genetic data suggested that 
the two refugial populations expanded around the Basin and 
re-connected at two zones. A broad contact zone at the south, 
which became an evolutionary 'melting pot' (Dufresnes et al 
2016; Qiao et al., 2018), incurred extensive admixture between the 
two refugial populations. The intermediate phenotypes of the 
southern populations are likely consequences of this admixture 
(Figure IB, Wilcoxon test on southern and western populations 
showed no significant). The second contact zone is located at 
the north-western corner of the Basin, where only limited gene 
exchange occurred. Concordantly, the morphological variations 
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also showed significant difference between the western 
and northern populations with limited intermediate forms 
(Wilcoxon test between northern and western populations 
were significant in both sexes: males, Wilcoxon W - 0.99, P « 
0.05; females, Wilcoxon W = 3.06, P « 0.005). 

It is probably difficult to completely reject one of the two 
alternative hypotheses, as they are not mutually exclusive and 
partially correlated. In the case of the green odorous frog, its 
evolutionary history combined with climatic factor clearly 
provides a better explanation of its morphological variation 
than either hypothesis itself. 

We study also rejected several other factors that often 
contribute to morphological variations. We tested the year 
of collection and locality altitude, and neither is correlated 
with the morphological traits that we examined (Table 3). 
Nevertheless, our small sample sizes may limit the detecting 
power of our analysis. Furthermore, we did not determine the 
ages of our samples. Since amphibians have indeterminate 
growth, the impacts of sample age need to be taken into account 
for formulating hypotheses in future studies. 

Odorrana margaretae has a ring-species alike genetic 
divergence pattern around the Sichuan Basin, which make 
it an excellent model system for studying speciation. This 
study revealed a compatible geographical clinal variation of 
its morphology. Future work on other phenotypic variations, 
particularly these may potential link to reproductive isolation, 
are essential to understand it species history and hybridization 
dynamics. 
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Appendix 


Table S1 Sample site and sample size information for specimens of the odorous frog Odorrana margaretae. 


Site Group Location nSamples nMales  nFemales Coordinates Elevations(m) Collected year 


N 


north Nanjiang, Guangwushan 22 19 3 106.40704E, 32.34705N 1412 2012 


EN 
= 
to 


north Chongqing, Wushan 16 4 109.97463E, 31.40330N 765 1956 


lon 
= 
© 
Re 


north Baokang, Wudaoxia 111.20576E, 31.71790N 519 2012 


co 
Ne) 
A 


east Sangzhi, Badagongshan 5 109.81263E, 29.67515N 1138 2011 


a 
o 
= 


10 south Hejiang, Zihuai 106.17358E, 28.37992N 763 2011 


o 
an 


12 south Bijie, Dafang 3 105.9991 1E, 27.54319N 852 2014 


— 


4 west Leshan, Mt. O'mei 11 


M 


4 103.38918E, 29.56418N 749 2012 


— 
M 


6 west Ya'an, Yanchang 16 9 103.07560E, 29.74135N 919 2012 


ja 
Un 
o 


8 west Dujiangyan, Hongkou 5 103.65208E, 31.12045N 1102 2012 


N 
u 
w 


0 west Pingwu, Laohegou 2 104.71170E, 32.46250N 1106 2018 
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BIO16 Precipitation of Wettest Quarter BIOIS8 Precipitation of Warmest Quarter 
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Figure $1 A: Patterns of geographical variation in precipitation related factors; B: Patterns of geographical variation in temperature 
related factors. For (A-B), mean values from every site were plotted in coordinates, and the grey scale of symbols is according to in- 
creasing size. 


Table S3 The mean PCI scores of PCA conducted on 21 morphological characters for individuals from each site. Males and females are 
shown separately. 


= 
£ 
a 


Site Group Female 


N 


north -0.71 -13.5 


A 


north 0.27 1.45 


lon 


north -11.11 


oo 


east -12.79 -11.46 


= 
o 


l 
N 
= 
w 


south 


— 
N 
o 
E 
B 


2.88 


I 
N 


— 
4A 


west 6.67 11.51 


= 
lon 


west 5.59 3.07 


=á 
oo 


west 27.04 


NX 
o 


west 5.56 3.63 


